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Summary 

The  exact  governing  equations  of  fluid  dynamics  are  too  computationally  expensive  to 
solve  on  a  computer  for  practical  applications.  They  will  remain  intractable  for  at  least 
the  next  50  years.  In  order  to  computationally  predict  the  performance  of  engineering 
applications  that  involve  fluids,  alternative  (and  computationally  tractable)  equations 
must  be  formulated.  These  equations  are  referred  to  as  turbulence  models. 

This  project  has  successfully  developed  an  entirely  new  theory  and  computational 
approach  to  modeling  turbulence.  The  idea  was  to  track  a  statistical  ensemble  of 
colliding  and  spinning  disks  (oriented  eddies)  in  colloidal  suspension  in  the  fluid.  This 
somewhat  usual  fluid  model  (akin  to  the  flow  ot  liquid  crystals)  was  chosen  because  we 
observed  the  remarkable  property  that  this  type  of  non-Newtonian  fluid  behaves  just  like 
a  turbulent  fluid  in  the  rapid  distortion  (large  mean  flow  gradients)  limit. 

The  project  was  successful  in  demonstrating  that  with  appropriate  extensions  this  model 
of  turbulence  could  be  applied  to  a  very  wide  variety  of  turbulent  flows  with  high 
predictive  accuracy.  A  number  of  different  limits  in  which  the  model  is  exact  were 
demonstrated.  We  were  also  successful  in  keeping  the  number  of  model  constants  very 
low.  In  addition,  the  theory  behind  the  model  was  vastly  expanded  when  we  were  able 
to  show  that  this  approach  is  equivalent  to  a  model  for  the  evolution  of  the  two-point 
velocity  correlation  tensor. 

The  prediction  (or  modeling)  of  turbulent  fluid  flow  is  arguably  one  of  the  greatest 
bottlenecks  in  the  Navy’s  ability  to  rapidly  design  innovative  devices  and  respond  to 
environmental  threats73.  While  this  research  does  not  address  a  specific  Navy  operational 
issue,  it  has  an  extremely  broad  reaching  impact  on  Navy  operations  in  general  and  will 
enhance  the  Navy’s  ability  to  successfully  execute  its  mission. 


Background 


One  of  the  greatest  bottlenecks  in  Engineering  Design  today  is  the  computational 
prediction  of  turbulent  fluids.  In  a  vast  variety  of  applications:  from  air  pollution,  to 
engine  efficiency  and  emissions,  to  global  climate  prediction,  to  submarine  performance, 
to  galactic  evolution,  turbulence  plays  a  critical,  if  not  the  most  important,  physical  role. 
Some  of  these  problems  will  never  be  computationally  tractable  without  a  turbulence 
mode!,  Those  that  do  become  tractable  in  the  coming  decades  (using  perhaps  large  eddy 
simulation)  will  then  represent  a  stunning  waste  of  resources  that  could  be  used  more 
productively  for  testing  and  optimizing  multiple  designs  or  adding  additional  physics. 

In  theory,  low  cost  turbulence  models  that  are  predictive  could  have  a  profound  effect  on 
how  computational  fluid  dynamics  (CDF)  is  used  in  the  design  process.  CFD  could  go 
from  being  a  qualitative  predictor  of  trends  to  a  pervasive  qualitative  design  tool 
comparable  to  Finite  Elements  and  computational  mechanics.  While  existing  engineering 
turbulence  models  do  not  currently  provide  predictive  accuracy,  it  is  demonstrated  in  this 
proposal  that  by  modeling  the  turbulence  structure  as  well  as  the  fluctuating  velocity 
magnitudes  the  Oriented-Eddy  Collision  (OEC)  model  is  able  predict  the  influence  of  the 
mean  flow  on  the  turbulence  exactly.  Since  this  is  the  dominant  physical  effect  on  the 
turbulence  evolution,  the  OEC  model  is  highly  predictive  (only  turbulence-turbulence 
interactions  require  modeling  and  this  process  is  relatively  well  understood).  The  price 
for  this  predictive  accuracy  is  an  order  of  magnitude  increase  in  the  cost  of  the  model 
compared  to  traditional  approaches.  This  puts  the  model  cost  at  a  few  times  the  cost  of 
the  mean  flow  calculation  (which  is  appropriate  given  the  critical  physics  contained 
within  the  turbulence)  and  many  orders  of  magnitude  less  expensive  than  large  eddy 
simulation  (LES)  or  direct  numerical  simulation  (DNS)  solutions.  This  is  a  region  of  the 
cost/performance  parameter  space  that  has  not  been  extensively  explored  previously  in 
the  context  of  turbulence  modeling  and  which  holds  great  promise. 

Significance  &  Objectives 

The  computational  resources  required  to  solve  the  governing  equations  of  fluid  dynamics, 
the  Navier-Stokes  equations,  increase  as  the  cube  of  the  Reynolds  number  (Rogallo  & 
Moin).  If  large  eddy  simulation  (LES)  models  are  eventually  successful  this  exponent 
could  be  reduced.  Still,  LES  of  a  passenger  jet  wing  is  not  expected  to  be  possible  on  the 
largest  supercomputer  until  around  2045  (Spalart).  This  assumes  current  increases  in 
computer  speed,  which  may  well  not  continue.  Most  environmental  problems  have  much 
larger  Reynolds  numbers  and  may  never  be  computationally  tractable  with  these 
approaches.  Even  problems  that  do  become  tractable  are  a  waste  of  resources  in  the 
context  of  Engineering  Design  where  Computational  Fluid  Dynamics  (CFD)  has  the 
potential  to  be  the  most  useful.  Computing  a  few  engineering  parameters  (such  as  wing 
drag)  with  LES  and  the  supercomputer  of  the  year  2045  is  akin  to  swatting  a  fly  with  a 
nuclear  bomb.  These  resources  could  be  used  much  more  effectively  to  optimize  the 
design  or  add  additional  physics.  Computers  are  having  a  profound  effect  on  engineering 
design,  but  progress  on  problems  involving  fluids  has  been,  and  will  remain,  stalled 
indefinitely  until  cost-effective  and  predictive  turbulence  models  are  developed. 
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The  purpose  of  this  project  was  to  develop  a  turbulence  model  that  included  more  physics 
(thereby  increasing  the  predictive  accuracy)  while  at  the  same  time  drastically  reducing 
the  number  of  empirically  determined  model  constants.  This  comes  at  a  significant 
computational  price,  but  not  a  price  that  is  out  of  line  with  the  importance  of  the  turbulent 
physics  that  the  model  is  required  to  represent.  It  is  a  price  which  is  many  orders  of 
magnitude  less  than  LES  or  DNS,  and  yet  still  provides  the  predictive  accuracy  necessary 
to  move  CFD  from  being  a  qualitative  predictor  of  trends,  to  a  pervasive  quantitative 
design  tool  comparable  in  utility  to  current  day  Finite  Elements  in  computational 
mechanics. 

Classical  Turbulence  Modeling 

The  coarse-grained  equations  for  incompressible  fluid  flow  are  given  by, 

=-P.:  +  -(»,«;  -«/“/),,  (la) 

(lb) 

These  look  almost  the  same  as  the  original  governing  Navier-Stokes  equations  except  for 
the  final  source  term  in  the  momentum  equation.  The  tensor  in  the  final  term  represents 
the  influence  of  the  unresolved  turbulent  fluctuations  on  the  coarse-grained  velocity 

evolution.  When  a  Reynolds  (ensemble)  average  is  used  Rt  =  u,u  -ufi,  is  referred  to  as 

the  Reynolds  stress  tensor.  In  large  eddy  simulation,  this  quantity  is  referred  to  as  the 
subgrid-scale  stress  tensor. 

The  Reynolds  stress  tensor  is  the  critical  unknown  that  results  from  coarse-graining  the 
Navier-Stokes  equations.  It  describes  how  turbulent  fluctuations  influence  the  mean  (or 
resolved/coarse)  flow  evolution.  Once  a  model  for  the  Reynolds  stress  tensor  is 
formulated  these  equations  are  much  easier  to  solve  computationally  than  the  original 
Navier-Stokes  equations.  Most  importantly,  when  mesh  adaptation  is  employed  the 
computational  cost  does  not  scale  with  Reynolds  number.  These  equations  look  like  the 
Reynolds  Averaged  Navier-Stokes  (RANS)  equations  formulated  over  a  century  ago 
(Reynolds),  However,  it  has  recently  been  recognized  that  these  equations  hold  for  any 
type  of  coarse  graining  (not  just  the  Reynolds  ensemble  average),  and  are  directly 
applicable  to  LES  (German o).  Since  the  fundamental  modeling  problem  and  equations 
are  the  same  for  both  RANS  and  LES,  the  model  discussed  in  this  proposal  is  equally 
applicable  (and  actually  very  well  suited  except  for  its  cost)  to  LES. 

In  the  context  of  turbulence  modeling  the  workhorse  of  both  commercial  and  in-house 
Computational  Fluid  Dynamics  (CFD)  codes  is  the  two-equation  models,  most  notably 
the  K-e  model.  These  models  involve  modeled  partial  differential  transport  equations  for 
two  scalar  turbulence  quantities  (most  often  the  turbulent  kinetic  energy,  K,  and 
dissipation  rate,  s).  There  is  also  a  hypothesized  algebraic  constitutive  equation  relating 
these  two  scalar  quantities  and  the  mean  (or  coarse)  flow  information  to  the  Reynolds 
stress  tensor.  In  many  respects,  the  predictive  performance  of  two-equation  models  is 
phenomenally  impressive  given  the  complexity  of  the  physical  process  that  they  are 
charged  with  representing.  They  work  quite  well  predicting  basic  turbulence  producing 
regions  such  as  simple  boundary  layers  and  free  shear  layers.  The  modifications  based 
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on  renormalization  group  theory  (Orzag  &  Yakhot)  and  elliptic  relaxation  (Durbin)  have 
even  expanded  the  predictive  scope  of  these  models.  Nevertheless,  it  is  well  understood 
at  this  time,  even  by  CFD  users,  that  the  physics  contained  in  these  models  is 
fundamentally  limited.  These  types  of  models  are  incapable  of  predicting  with  real 
predictive  accuracy  the  wide  variety  of  turbulent  conditions  found  in  real  engineering 
design  applications. 


it  was  expected  that  Reynolds  Stress  Transport  (RST)  models  would  perform 
considerably  better.  These  models  do  not  use  a  hypothesized  constitutive  equation  for  the 
Reynolds  stress  tensor.  Instead,  they  start  with  the  exact  (but  unclosed)  transport 
equation  for  the  Reynolds  (or  subgrid  scale)  stress  tensor, 


^IJJ  ^ik  Ji  i 
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(2) 


The  last  three  source  terms  (the  transport  term,  the  pressure  term,  and  the  dissipation 
term)  require  modeling.  The  lack  of  an  ad  hoc  constitutive  equation  is  a  very  important 
step,  but  just  as  important  is  the  fact  that  the  influence  of  the  mean  flow  on  the  turbulence 
evolution  (the  second  term  on  the  right  hand  side  of  the  transport  equation)  is  closed  and 
does  not  require  a  model.  Mean  flow  gradients  are  the  primary  mechanism  for  turbulence 
generation  and  this  equation  suggests  that  this  mechanism  is  exactly  captured. 


Unfortunately,  the  predictive  promise  of  RST  models  has  not  been  realized  and  while 
RST  models  are  available  in  most  CFD  codes  they  are  rarely  used.  In  practice,  their 
predictive  performance  is  rarely  better  than  the  more  computationally  robust  two- 
equation  models.  The  problem  lies  in  the  fact  that  the  pressure  term  also  carries 
information  about  the  mean  flow  gradients  which  is  just  as  important  as  the  explicit 
production  term  but  which  is  extremely  difficult  to  model.  The  effect  of  the  mean  flow 
gradients  is  therefore  not  exactly  captured.  The  modeled  pressure  term  must  capture 
those  effects.  It  has  even  been  shown  that  modeling  this  term  is  fundamentally 
impossible  in  certain  circumstances  using  only  the  information  provided  by  the  RST 
models  (Reynolds,  Speziale).  In  addition  to  this  fundamental  shortcoming  RST  models 
have  another  significant  source  of  error  and  uncertainly.  This  is  the  ad  hoc  transport 
equation  (usually  an  e  equation)  that  is  required  to  close  the  system.  Other  issues  exist 
but  these  two  effects  (the  pressure-term  and  additional  transport  equation)  represent  the 
majority  of  empirical  model  constants  and  predictive  uncertainty  in  RST  models. 


The  Problem  with  Classical  Models 

Some  typical  results  for  a  RST  modei  are  shown  in  Figure  1 ,  In  each  of  these  situations, 
homogeneous  isotropic  turbulence  is  subjected  to  a  very  simple  arrangement  of  mean 
flow  gradients.  Depending  on  the  particular  arrangement  of  mean  flow  strains,  certain 
components  of  the  Reynolds  stress  tensor  are  amplified  with  time  and  others  decrease.  In 
Figure  1,  the  circles  are  direct  numerical  simulation  (DNS)  data  and  the  lines  are  the 
corresponding  RST  model  predictions.  In  some  cases,  the  model  predictions  are  quite 
accurate;  however,  in  others  the  agreement  is  poor.  The  level  of  predictive  performance 
can  be  characterized  by  the  dimensionless  velocity  gradient,  S*  =(ut  ^  t)'12 K / e  .  This 
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(a)  (b)  (c) 


(d)  (e)  (0 


Figure  l:  State-of-the-art  RST  model  predictions  of  turbulence  subjected  to  simple  mean  flows. 
Dimensionless  normalized  shear  stress  as  a  function  of  dimensionless  time.  Symbols  are 
experimental  or  DNS  data  and  lines  are  RST  model  predictions,  (a)  plain  strain,  S*  =  0.5  (b) 
plain  strain,  S*  =  4  (c)  axisymmetric  contraction,  S*  =  4.08  (d)  axisymmetric  expansion,  S*  = 
55.8  (e)  homogeneous  shear,  S*  =  4.71,  (f)  homogeneous  shear,  S*  =  30.9.  Data  of  cases  (a-d) 
are  from  Lee  &  Reynolds  and  (e-f)  from  Matsumoto. 

is  a  ratio  of  the  turbulence  response  time  to  the  mean  flow  time  scale,  If  S*  is  small  the 
turbulence  response  time  is  much  less  than  the  time  over  which  the  mean  flow  changes 
appreciably.  In  this  case,  the  turbulence  is  always  in  quasi-equilibrium  with  the  mean 
flow  (since  it  can  respond  much  more  quickly  than  the  mean  flow  is  changing)  and  the 
RS I  model  performs  well.  However,  when  the  turbulence  is  subjected  to  rapid  changes 
the  model  performs  poorly.  RST  models  are  fundamentally  quasi -equilibrium  models. 

There  are  many  flow  situations  such  as  turbulent  boundary  layers  and  shear  layers  where 
the  flow  is  in  quasi-equi librium.  However,  in  engineering  applications  there  are  usually 
critical  flow  regions  where  quasi-equilibrium  is  not  a  good  approximation.  These  include 
separating  flows,  transition,  reattaching  flows,  fast  moving  boundaries  (such  as  turbine 
blades)  and  flows  under  the  influence  of  strong  body  forces  such  as  gravity  or 
electromagnetic  forces.  While  the  non-equilibrium  regimes  occupy  less  volume  in  a 
typical  C’FD  calculation  that  the  quasi-equilibrium  regimes  they  often  dominate  the 
critical  physics  and  frequently  dictate  the  overall  predictive  performance  of  the 
turbulence  model.  Figure  2a  shows  a  pictorial  representation  of  parameter  space. 
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Accurate  RST 
or  k/s  models 


Poor  model  performance 


Accurate  OEC 


(a) 


(b> 


Figure  2:  Pictorial  representation  of  state-space  and  the  regions  of  applicability  of  (a)  traditional 
engineering  turbulence  models  and  (b)  the  eddy  interaction  model.  Traditional  models  perform 
well  close  to  the  origin  where  the  turbulence  is  in  quasi-equilibrium  with  the  mean  flow  and/or 
external  forces.  The  OEC  model  increases  the  range  of  predictive  performance  significantly  by 
exactly  capturing  mean  flow  (or  external  force)  effects  that  dominate  when  the  flow  is  not  in 
equilibrium. 


Traditional  turbulence  models  are  predictive  only  near  the  origin  of  this  diagram  where 
the  flow  is  close  to  some  equilibrium  state. 

What  is  most  interesting  is  that  the  highly  non-equilibrium  regime  should  be  the  easier 
regime  to  model.  A  high  degree  of  nonlinearity  implies  a  strong  separation  of  time 
scales,  which  in  turn  implies  that  course  graining  could  be  accomplished  very  accurately. 
While  turbulence  does  not  contain  a  separation  of  length  scales  (the  large  turbulent  length 
scales  are  always  the  same  order  of  magnitude  as  the  mean  flow  length  scales),  it 
sometimes  does  contain  this  time  scale  separation  that  can  be  advantageous.  The  OEC 
model  described  below  is  exact  in  the  limit  of  very  strong  mean  flow  gradients  (the  non- 
equilibrium  limit).  As  shown  in  Figure  2b,  this  implies  that  very  far  from  the  origin  of 
the  diagram  the  model  is  exact.  This  type  of  model  is  inherently  much  more  predictive. 

The  influence  of  the  mean  flow  on  the  turbulence  (or  Reynolds  stress  tensor)  evolution  is 
a  linear  process  (only  turbulence-turbulence  interactions  are  nonlinear).  Therefore,  in 
theory,  it  should  be  possible  to  capture  this  effect  exactly,  and  this  is  indeed  the  case. 
The  key  insight  that  is  necessary  to  do  this  is  to  recognize  that  although  the  Reynolds 
stress  tensor  contains  sufficient  information  to  describe  the  turbulence  effects  on  the 
mean  flow  it  does  not  contain  enough  information  to  correctly  represent  the  influence  of 
the  mean  flow  on  the  turbulence1  ’.  To  capture  the  influence  of  the  mean  flow  exactly  the 
model  must  also  capture  the  effect  of  the  mean  flow  gradients  on  the  turbulence  structure. 


The  Solution 

In  this  work,  we  use  two-point  correlations  as  a  measure  of  turbulence  structure.  If  two 
nearby  points  in  a  physical  domain  have  a  very  high  correlation  between  their  velocities. 
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this  implies  that  turbulence  structures  tend  to  be  aligned  with  the  line  joining  those  two 
points.  If  two  points  have  a  low  or  negative  velocity  correlation,  it  implies  that  the 
turbulence  structure  does  not  span  the  distance  between  those  two  points. 


One  major  result  of  this  project  has  been  the  discovery  that  the  OEC  model  is  related  to 
the  exact  two-point  correlation  equation. 


C>j.i  (*> r)  =  ~ui,k(x)Cki  ~ u/.k (*)Q  "  ut +  (u* (*) " «* 


«it»  \ 

I  dpu) 

Spu)  \ 

-®“k+v| 

f  s% 

9  ^ y 

\ 

r\  / 

drk 

\  * 

*i  J 

Sr,  T  V  | 

yS**31*  ’ 

^  J 

where  r  =  x-xis  the  distance  between  two  points  in  the  domain.  The  pressure  velocity 
correlation  is  actually  closed  in  this  equation.  It  does  not  require  a  model  and  can  be 
obtained  from  the  equation, 

2  (4) 


Only  the  triple  velocity  correlation,  Tl{Jk) ,  representing  turbulence-turbulence  interactions, 

requires  a  model  to  close  Eqn  (3).  The  Reynolds  stress  tensor  (needed  to  predict  the 
mean  flow)  is  given  by  RtJ  (x)  =  C(J  (x,  0) . 


The  reason  this  equation  has  not  been  used  as  a  starting  point  for  practical  turbulence 
models  in  the  past  is  due  to  its  high  dimension.  A  very  coarse  discretization  of  r,  with  a 
10x10x10,  mesh  implies  the  need  to  solve  1000  transport  equations  for  each  point  in 
physical  space.  This  cost  is  too  high.  A  key  cost  savings  comes  from  recognizing  that 
only  very  basic  information  about  how  the  correlations  behave  in  r  is  necessary  to 
capture  the  Reynolds  stress  evolution  (at  r  =  0)  correctly.  Our  real  goal  is  to  model  the 
Reynolds  stress  evolution  accurately,  not  the  two-point  correlation,  We  therefore 
parameterize  the  shape  of  the  correlation  with  the  vector  parameter,  k.  and  assume  that 
the  sum  of  a  few  parameterized  correlations  is  sufficient  to  represent  a  general  correlation 
function, 

e,(i.r)-£  4,00«'M  (5) 

Experimental  data  suggests  that  an  exponential  is  a  good  parameterization  for  the 
correlation  function  (Compte-Bellot  &  S.  Corrsin).  We  have  briefly  explored  the  use  of 
other  functions  and  this  aspect  of  the  model  will  be  further  explored  in  the  proposed 
work.  The  model  might  be  considered  to  be  an  assumed  shape  two-point  correlation 
model.  It  gives  roughly  the  right  balance  between  physical  accuracy  and  computational 
cost. 


Inserting  the  parameterized  two-point  correlation  (Eqn.  5)  into  the  exact  two-point 
correlation  equation  (Eqn.  3),  and  making  a  similar  decomposition  for  the  velocity- 
pressure  correlation,  results  in  the  following  equations. 


Kjt  -~kkuk,i  Lk,+ml 


(6b) 
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where  ±  =  (Kk2)V2  is  the  inverse  turbulence  timescale.  The  terms  a„,  D  ,  m,  represent 
the  nonlinear  turbulence-turbulence  interactions  and  must  be  modeled. 

The  models  for  both  DtJ  and  m,  each  currently  involve  one  empirical  model  parameter. 

However,  as  a  part  of  the  proposed  work  plan  below  we  have  definitive  approaches  that 
could  reduce  the  number  of  empirical  model  constants  to  zero  for  these  two  terms.  The 
inverse  timescale  constant  aH  was  set  via  DNS  data  for  isotropic  decay.  It  defines  the 
cross  over  between  high  and  low  Reynolds  numbers  and  has  very  little  influence  on  most 
predictions  (which  are  at  high  Reynolds  numbers).  A  good  parameterization  of  the  two- 
point  correlation  requires  a  number  of  these  equation  systems  to  be  solved  (and  then 
summed,  Eqn  5).  The  exact  number  of  systems  is  a  question  to  be  explored  in  more 
detail  in  this  project  but  preliminary  results  are  presented  below  in  the  cost  section  of  the 
proposal. 

It  has  been  shown  analytically  and  is  demonstrated  in  some  of  the  figures  below  that  this 
model  gives  exact  results  in  the  limit  of  large  mean  velocity  gradients.  In  addition,  the 
model  captures  both  very  low  Reynolds  number  and  very  high  Reynolds  number 
isotropic  decay  exactly.  Isotropic  decay  depends  on  nonlinear  turbulence-turbulence 
interactions  (the  energy  cascade).  While  this  process  must  be  modeled,  it  is  notable 
there  are  no  empirical  constants  associated  with  the  cascade  model  in  the  eddy  collision 
approach. 

The  OEC  model  is  almost  unique  in  its  ability  to  capture  mean  flow  effects  on  turbulence 
exactly.  The  other  two  examples  of  models  with  this  ability  are  based  on  the  Fourier 
transform  of  the  Navier-Stokes  equations  (Reynolds  &  Kassinos,  Van  Slooten  &  Pope). 
They  appear  to  be  strictly  limited  to  homogenous  turbulence  and  periodic  domains.  The 
two-point  equation  (Eqn.  3)  and  the  resulting  OEC  model  are  fully  general  and  have  no 
such  restrictions. 


Detailed  Equations 


Two  main  equations  are  used  to  represent  the  oriented  collision  model.  The  first  one 
represents/?^,  which  is  the  Reynolds  stress  (average  velocity  fluctuations)  for  one 

orientation  kt  (see  equation  7b  below).  The  orientation  vector, A:,,  has  units  of  1/length 
and  captures  the  eddy  size  and  orientation. 

K,  =  K  [«a  +  2<>  7T  -  Sa  )1  +  K  r +  2 u,  k  (M  -  %  )1 

_  (7a) 

~(aLvk2+aH^)RtJ-aH(^-^nl+(Rlj^+^)ml+V(y  +  vT)VRIJ 

where  (77)  =  (77)“  K'  ifiPj  KKf)^-  The  tota'  Reynolds  stress  defined  as  R;j  is  the 


averaged  sum  of  the  individual  R,r  meaning  R.  =jjI.Rlr .  Equation  (7a)  has  seven  grouped 
terms.  The  mean  flow  gradients  and  system  rotation  is  accounted  for  by  u‘k  =ul  k  +  etk  Qt , 
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with  Oj  being  the  rotation  vector  for  a  non-inertial  frame.  The  dissipative  behavior  of  the 
model  is  captured  by  [a^k1  +aH  -d-)/?u  andorw^Ji)y  is  the  retum-to-isotropy  model 
discussed  in  section  3.3  below.  The  factor is  the  timescale  used  to  model  the 
dissipation.  -^nt  is  the  rotation  term.  The  six  term  ( Rt)  p  +  kh  arises  from  the  need 

to  maintain  orthogonality  (Ryk,  =0)  between  the  orientations  and  the  RtJ  ( m,  is  the  k- 

return  model).  Incompressibility  requires/?,^,  =0.  The  final  term  V  (V  +  vr) VR:j  models 
the  diffusive  action  of  the  Reynolds  stresses. 


The  second  equation  represents  the  orientation  k,  with  its  time-derivative  defined  as: 

K>  =-**“*.,  ~Ya,yk2  +a“  ~^k‘  +in'  +(iH  +  V(u  +  vr)Vk, 


(7b) 


The  above  equation  contains  six  grouped  terms.  The  first  terms  captures  the  mean 
gradient  effects  (shear).  This  term  is  the  equation  for  passive  disks.  Just  as  in  (2.3.1),  the 

second  term  captures  the  dissipation;  l  takes  on  the  value  3  or  5  respectively  for  the  k 
or  k  low  wave  number.  The  third  term  n,  present  in  models  the  secondary  rotation 

effects  and  m,  is  the  return  model  for  the  orientations.  The  last  term, 
V(v  +  vy)V&,  accounts  for  the  diffusive  action  of  the  orientation  vectors/:, . 


In  addition. 

K.> = Z  At  -  2“,‘*  ]  -  «^Z  -  ^ZGrH 

(8) 

Hence, 

(9) 

Ability  to  Represent  Two-point  Correlations 

The  unknowns  in  the  oriented  eddy  model  are  closely  related  to  the  two-point 
correlations.  In  this  section,  we  take  a  brief  look  at  this  relationship.  Assume 

RtJ(xJ)  ~  'ERl)F{k  -r)  where  F(tj)  is  a  simple  function  of  r  ,  the  distance  between  two 

points.  Considering  the  specific  case  where  F(k,r )  =  e  ,  then  Rtj  (x,r)  =  £ • 
When  looking  at  the  two  point  correlation  in  the  x-direction  for  example,  we  get 

Similarly  for  R22  and  R33: 
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Hence,  we  obtained  the  contour  plots  shown  below  that  look  very  similar  to  the  two-point 
correlations  found  from  DNS  data. 


Figure  3.  Rn.  R22  and  R33  as  seen  from  the  redirection. 


DNS  two-point  correlation  data  corresponding  to  the  first  two  figures  above  (the  model) 
is  shown  below  in  Figure  4,  The  shapes  are  very  similar.  The  mesh  size  used  in  the 
DNS  simulation  was  768  by  768  by  1536  cells,  with  a  domain  size  of  56.54  by  56.54  by 
1 13.09  centimeters. 
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Figure  4.  A  planar  slice  of  a  three  dimensional  Rn  and  R22  two-point  correlation  in  the  X-Y 
plane  about  Z=0. 


A  major  break  through  of  the  current  modeling  approach  is  the  ability  to  model  the 
‘blocking’  affect  of  walls  exactly.  After  a  vertical  wall  is  inserted  (  at  y  =  -2  in  the 
figures  above)  the  two-point  correlations  are  dramatically  altered.  This  is  shown  in 
figure  5. 
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Figure  5.  A  planar  slice  of  Rt  1  and  R22  two-point  correlations  in  the  X-Y  plane  about  Z=0  after  a 
wall  is  inserted  at  y=-2. 


The  essence  of  the  exact  wall  representation  is  shown  in  figure  6  for  the  R22  correlation 
This  correlation  is  chosen  to  illustrate  the  idea  because  it  is  the  most  dramatically  altered 
correlation.  In  the  figure  below,  the  red  is  DNS  data  and  the  blue  is  the  model.  This 
figure  is  the  R22  correlation  along  a  line  perpendicular  to  the  wall.  Only  results  to  the 
right  of  y— -2  matter.  To  the  left  of  y=-2  there  is  now  a  solid  wall.  The  model  uses 
image  eddies  to  exactly  capture  the  influence  of  the  wall  on  the  two-point  correlations. 


I-  igure  6.  A  demonstration  of  how  wall  effects  can  be  exactly  captured  by  the  eddy  OEC  model. 
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NUMERICAL  RESULTS 


Below  is  the  table  that  summarizes  the  different  sections  and  results  of  the  current 


Table  1.  Research  Summary 


1.  Isotropic  Decay 

In  general,  when  the  properties  of  a  material  are  the  same  in  all  directions,  the  material  is 
said  to  be  isotropic.  In  the  case  of  turbulence,  if  the  fluctuations  are  independent  of 
direction,  the  turbulence  is  isotropic.  When  the  fluctuations  do  not  have  any  directional 
preference,  then  the  off-diagonal  components  of  R(J  vanish,  and  Rt]  =  R22=R,3. 

Mathematically,  this  corresponds  to  Ry=|K  . 

In  this  work,  it  is  necessary  to  define  isotropy  for  the  orientations  as  well.  For  isotropy, 
all  orientation  vectors  have  the  same  magnitude  and  arc  uniformly  distributed  on  the 
sphere. 
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Von  Karman  &  Howarth  first  suggested  in  1938  that  the  decaying  turbulence  should  have 
a  power  law  behavior  of  the  form: 


K  =  K„ 


l+3L 


nK 


0  / 


(10) 


where  K0  is  the  initial  turbulent  kinetic  energy  and  e0  represents  the  initial  dissipation, 
and  n  is  the  decay  exponent.  While  all  researchers  agree  on  the  power  law  form,  there  is 
less  agreement  on  what  the  value  for  n  should  be.  However,  most  investigators  agree  that 

the  exponent  n  is  highly  dependent  on  the  low  wavenumber  k  of  the  energy  spectrum 

(Saffman).  In  the  case  where  the  low  wavenumber  portion  of  the  spectrum  goes  as  A2 ,  n 
corresponds  to  3/2  at  low  Reynolds  number  and  6/5  at  high  Reynolds  number.  On  the 

other  hand,  when  the  low  wavenumber  portion  of  the  spectrum  goes  as  A4 ,  n  corresponds 
to  5/2  for  low  Reynolds  number  and  10/7  for  high  Reynolds  number. 


We  will  attempt  to  obtain  all  these  limits  with  the  OEC  model.  For  isotropic  decaying 
turbulence,  the  dissipation  s  is, 
dK 

s=-j  <n> 

Substituting,  equation  (3. 1.1,1)  to  (3. 1.1.2)  above,  we  obtain: 


£  =Sn 


1+  £f>t 


\ 


-n-I 


nK 


0  J 


02) 


In  this  section,  our  model  attempts  to  capture  the  evolution  of  n  as  a  function  of  the 


£ 

turbulent  Reynolds  number  ( Rer  —  — ,  with  v  being  the  fluid  kinematic  viscosity)  for 

low  wavenumber  (not  orientations)  of  A2  and  A4 .  Figure  4  below  summarizes  the  results 
obtained  when  the  low  wavenumber  behavior  of  the  spectrum  is  A2. 
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RCy 


A  _ 

Figure  7;  Power-law  exponent  as  a  function  of  the  turbulent  Reynolds  number  for  a  k 2  low 
wavenumber  spectrum. 


The  thick  maroon  and  dark  green  lines  represents  the  mode!  predictions  (a,  =6,15,30 
for  aH  =  1 ).  For  our  purpose,  we  determined  that  the  ratio  —  =  1 5  (maroon  curve)  best 

<*H 

matched  the  DNS  simulations  of  Chasnov,  Mansour  &  Wray.  The  upper  and  lower  purple 
dashed  lines  included  in  the  figure  are  the  low  and  high  Reynolds  bounds  on  n.  Notice 
that  the  model  obtains  these  limits  independent  of  al  .  Also  on  Figure  7  are  shown  the 
exponent  values  for  the  DNS  of  de  Bruyn  Kops. 

When  the  low  wavenumber  behavior  of  the  spectrum  goes  as£4 ,  we  obtained  the  results 
shown  in  Figure  8  below.  Again,  the  horizontal  green  dash  lines  represent  the  upper  (5/2) 

A 

and  lower  (10/7)  limits  of  the  exponent  for  &  spectrum.  The  thick  purple  and  blue  lines 

X 

are  the  model  predictions  for— =  10,25,50 .  In  addition  to  these  curves,  there  are  four 

XH 

I283  DNS  simulations  by  Yu  et  al.  and  four  256}  DNS  simulations  by  Mansour  &  Wray. 

X 

For  the  same  reason  mentioned  above,  we  determined  that  —  =  1 5  (not  shown)  is  an 

XH 

adequate  compromise.  Note  that  this  ratio  is  similar  to  the  one  determined  above  for  k 2 . 
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Figure  8.  Power  law  exponent  as  a  function  of  Reynolds  number  for  a  kA 

spectrum. 


low  wavenumber 


Kinetic  Energy 

In  this  section,  we  focus  our  efforts  on  predicting  the  decay  of  kinetic  energy  in  isotropic 
flows  (other  than  just  the  exponent).  This  is  essentially  a  posteriori  test  of  the  chosen 
A. 

value  — ^  =  15*  We  test  the  model  against  numerous  published  data:  some  experimental, 

Ah 

some  LES  and  other  DNS,  In  determining  the  kinetic  energy,  the  equations  used  in  our 
model  predictions  originate  from  equation  (7a)  and  (7b)  above  with  the  particularity  that 
the  flow  is  isotropic.  Hence,  there  is  no  need  to  include  the  return-to-isotropy  ( Dtl  =  0 )  as 

well  as  the  diffusion  terms:  V(i/  +  vy7,)V^y  =  0,  V(i/  +  vr)V&(  =0  .  Thus,  in  cases  where 
no  rotation  is  present, 

k.'=-{l5vk'+-t)k  (13a) 

k,,-~(L5vlc,+-j;)k,  (l3b) 

where  /  =  3  for  k~  and  /  =  5  for  kA  low  wavenumber  spectra. 
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Data  and  model  predictions  are  shown  below  for  low  and  intermediate  turbulent 
Reynolds  numbers,  In  addition,  we  state  all  initial  conditions  in  Table  2  below: 


Wigeland  &  Nagib 

(exp.  Data) 

Man 

Carol 

Spe: 

(D1 

sour* 

3on  & 
dale 

MS) 

Jsacquin 

{exp.  Data) 

de 

Bruyn 

Kops 

(DNS) 

Squires 

(LES) 

14.85 

2.96 

2.77 

0.93 

0.95 

11.73 

16.43 

30.93 

0.782 

1.27 

1.35 

K(mW) 

0.098 

0.045 

0.029 

0.964 

0.977 

0.15 

0.264 

0.462 

0.087 

0.265 

0.298 

v{mVs) 

1.8 

1.8 

1.8 

3.67 

1.49 

1,51 

1.51 

1.51 

1.49 

8.6 

8.0 

e-5 

e-5 

e-5 

e-2 

'  e-2 

e-5 

e-5 

e-5 

e-5 

e-5 

e-5 

RCr 

36 

38 

17 

27.2 

67.1 

127 

281 

457 

655 

643 

764 

Table2,  Initial  conditions 


In  Figure  9,  the  kinetic  energy  is  represented  versus  time.  The  asterisks,  the  triangles  and  the  stars 
correspond  to  the  experimental  data  with  corresponding  Rex=36,  38  and  1 7  while  the  dashed  lines 
correspond  to  our  simulations. 
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Figure  9.  Wigeland  and  Nagib's  decaying  kinetic  energy. 

In  Figure  1 0,  the  kinetic  energy  versus  time  is  shown.  The  orange  dots  correspond  to  the 
experimental  data  for  Rex=27.24  and  the  purple  ones  are  for  Rer=67.1.  The  solid  lines 
correspond  to  our  simulations.  Clearly,  the  OEC  model  shows  good  agreement  with  the 
DNS  data  of  Mansour,  Cambon  and  Speziale. 
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Figure  10.  Mansour,  Cambon  and  Speziale’s  decaying  kinetic  energy. 


Figure  1 1  shows  the  kinetic  energy  versus  time.  The  asterisks  correspond  to  the 
experimental  data  and  the  dashed  lines  correspond  to  the  simulations. 
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Figure  1 1.  Jacquin’s  decaying  kinetic  energy. 

Figure  12  shows  the  kinetic  energy  versus  time.  The  red  asterisks  correspond  to  DNS  data 
of  de  Bruyn  Kops  for  Ref=655  and  the  dashed  lines  correspond  to  the  oriented  eddy 
model  simulations. 


Figure  12.  de  Bruyn  Kops’s  decaying  kinetic  energy. 
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figure  13  shows  the  kinetic  energy  versus  time.  The  green  asterisks  correspond  to  the  k2 
experimental  data  with  Rer  =  643  .  The  red  asterisks  represent  k*  data  with  Rer  =  764 . 
The  blue  and  pink  lines  correspond  to  our  simulations. 


Figure  13.  Squires’  decaying  kinetic  energy  for  both  k1  and  &4 . 


Based  on  the  data  presented  above,  it  is  concluded  that  the  OEC  model  performs  well  in 
predicting  the  decaying  kinetic  energy  for  simple  (homogeneous,  isotropic  and 
irrotational)  turbulent  flows. 


Rotating  Decaying  Grid-Turbulence 


To  measure  the  degree  of  rotation  present  in  the  flow,  we  used 
number  defined  as  follow; 


a  turbulent  Ross  by 
(14) 
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Large  Rq7  means  no  rotation,  whereas  Ror  <  i  implies  a  flow  dominated  by  rotation. 
With  rotation  present,  the  model  equations  become: 

*u=-H15v*2 +-£)*, -£«,  (15) 

Three  models  for  the  rotation  term  were  tested: 

A  Ml 

n  = - - — ! - -k 

(C^+c^Kn'D 

or 


«.*=■ 


(16) 


or 


1  (c1*‘i:+c1(n,)J)  f 

where  Q,  =  £<JkUkJ  +  Q,fnm" .  In  earlier  work  Chartrand  briefly  looked  at  the  first  two 

models,  nA  and  n,B  .  However,  after  extensively  studying  the  performance  of  each  of 

these  models  and  comparing  them  to  multiple  DNS  results,  we  came  to  the  conclusion 
that  the  above  two  terms  each  only  captures  a  different  aspect  of  the  rotation.  Hence,  the 
third  model  was  developed. 


With  each  model,  come  two  constants  C,  and  C2  that  are  used  to  tune  the  model 
behavior.  That  is,  C,  and  C2  are  both  model-dependent.  From  equations  (3. 1.3. 3)  above, 
it  is  clear  that  C2  affects  simulations  at  large  rotation  rates  while  C,  acts  at  small  rotation 
rates.  We  used  this  concept  in  determining  the  values  for  both  C,  andC2.  Table  3  below 
summarizes  the  values: 


Model 

Formula 

c, 

c2 

k 

k 

8 

0.25 

Q 

,  o' 

(c,*!ir+ci(f i*)  >  1 

20 

Va 

Smooth  k 

em2  , 

(C1*,*+CI(  n*)1)  ' 

20 

% 

Table  3:  Rotation-models  along  with  their  respective  tuning  constants C,  and  C, 

Next,  we  compared  the  performance  of  each  model  for  three  sets  of  data:  Jacquin, 
Mansour  Cambon  &  Speziale,  and  Blaisdell.  The  k-smooth  model  outperforms  the  other 
two.  The  initial  conditions  are  shown  in  Tables  4.  5  and  6  below: 
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Mansour,  Cambon  &  Speziale 

Jacquin 

Blaisdell 

0.93 

0.95 

11.73 

16.43 

30.93 

1.78 

K(m2/s0 

0.964 

0.977 

0.153 

0.288 

0.444 

1 

v(m2/s) 

3.67e-2 

1.49e-2 

1.5  le-5 

!.51e-5 

1.5  le-5 

4.4Ie-2 

Rer 

27.2 

67.1 

127 

281 

457 

12.75 

Rot 

0.37 

0.037 

0.24 

0.1 

1.22 

0.91 

1. 10 

— 

s 

— 

— 

— 

— 

— 

— 

— 

3 

Table  4:  Initial  conditions  of  Mansour,  Cambon  &  Speziale,  Jacquin  and  Blaisdell, 


Wigeland 

&  Nagib 

£<m2/sJ) 

14.67 

14.94 

3.49 

3.36 

3.36 

22.26 

K(m2/s2) 

0.0975 

0.105 

0.0462 

0.051 

0.033 

0.096 

v(m2/s) 

1.8e-5 

1.8e-5 

1.8e-5 

!.8e-5 

1.8e-5 

1.8e-5 

R©t 

36 

41 

34 

43 

18 

23 

Rot 

7.52 

1.78 

3.77 

0.82 

5.09 

2.9 

Table  S:  Wigeland  &Nagib’s  initial  conditions. 


s 

himomura 

de  Bruyn  Kops 

Veeravalli 

e(m  V) 

0.024 

0.025 

0.028 

0.0992 

7.96 

8.13 

0.098 

0.2619 

0.5638 

5.888e-2 

0.17 

0.202 

v(m2/s) 

r  8.0e-3 

8.0e-3 

8.0e-3 

1.4854e-5 

1 ,6e-5 

1.6e-5 

Rer 

50 

343 

1419 

2353 

227 

313 

Rot 

N/A 

0.095 

0.017 

0.006 

0.5 

0.32 

Table  6:  Initial  conditions  of  Shimomura,  de  Bruyn  Kops  and  Veeravalli. 


In  figure  14  below,  the  performance  of  each  model  is  analyzed  using  the  DNS  data  from 
Jacquin.  Our  simulation  matched  the  dimensionless  initial  conditions  of  Jacquin  as 
represented  in  Table  4.  The  crosses,  stars  and  dots  represent  the  experimental  data.  The 
solid  lines  represent  co,  the  dashed  lines  k  and  the  dotted  lines  k-smooth. 
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Figure  14:  Performance  comparison  of  k,  a>  and  k-smooth  rotation  terms. 

Looking  at  the  graph  above,  it  is  concluded  that  ail  three  rotation  models  performed 
equally  in  this  case,  due  to  the  somewhat  identical  turbulent  Ross  by  numbers  (1 .22,  0.91 
and  1.10). 
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In  Figure  15,  the  dimensionless  initial  conditions  of  Mansour,  Cambon  &  Speziale  were 
matched  for  Rei=27.24. 


Figure  1 5.  Performance  comparison  of  k,  o)  and  k-smooth  rotation  terms  based  on  Mansour, 
Cambon  and  Speziale  experimental  data,  a)  Ro=0.37.  b)  Ro=0,037. 
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In  Figure  16,  the  dimensionless  initial  conditions  of  Mansour,  Cambon  &  Speziale  were 
matched  for  Rej=67,l. 


Figure  16:  Performance  comparison  of  k,  o)  and  k-smooih  rotation  terms  based  on  Mansour, 
Cambon  and  Speziale’s  experimental  data,  a)  Ro=0.24.  b)  Ro=0.l . 
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We  also  evaluated  all  three  rotational  models  for  homogeneous  flows;  specifically,  using 
data  from  Blaisdell’s  elliptical  flow  as  shown  below  in  Figure  17: 


Figure  17:  Performance  comparison  of  k,  co  and  k-smooth  rotation  terms  for  Blaisdell, 

(homogeneous  shear  flow) 


From  the  graphs  above,  the  k -smooth  model  is  always  consistently  between  the  k  and  the 
( o-models .  And  sometimes  the  difference  is  so  subtle  that  it  is  almost  negligible.  In  the 
Blaisdell  case  however,  the  w  co-model  performs  very  poorly.  Hence,  it  was  decided  that 
the  k-smooth  model  performs  the  best.  So,  the  OEC  model  was  tested  against  other 
published  data  such  as  Wigeland  &Nagib,  Jacquin,  Shimomura,  de  Bruyn  Kops, 
Veeravalli  and  Mansour,  Cambon  and  Speziale.  (The  initial  conditions  are  presented  in 
Tables  5  and  6  above).  The  non-rotating  initial  conditions  for  Wigeland  &  Nagib  as  well 
as  de  Bruyn  Kops  were  already  given  above  in  the  prior  section  on  decaying  turbulence. 

In  Figure  1 8  below,  the  asterisks,  triangles  and  stars  represent  the  experimental  data  of 
Wigeland  &  Nagib  for  low  Reynolds  number,  while  the  dotted  lines  represent  the 
predictions  for  the  collision  model. 
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Figure  18:  Rotating  isotropic  decay  of  Wigeland  &  Nagib  using  the  rotation  model  n‘  . 

Turbulent  kinetic  energy  versus  time. 


In  Figure  19,  the  asterisks,  crosses  and  squares  represent  the  experimental  data  of 
Jacquin,  while  the  dashed  lines  represent  the  predictions  of  the  collision  model  with 
n-  for  rotation  model.  The  numbers  in  parenthesis  140,  310  and  500  correspond  to  the 
turbulent  Reynolds  number. 


Figure  19:  Rotating  Isotropic  decay  of  Jacquin.  Turbulent  kinetic  energy  versus  time. 
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In  Figure  20,  the  asterisks  represent  the  experimental  data  of  Shimomura  (for  both 
irrotational  and  rotational  cases),  while  the  solid  lines  represent  the  predictions  from  our 
collision  model.  As  summarized  in  Table  6  above,  the  turbulent  Reynolds  numbers 

correspond  respectively  to  50,  343  and  1419.  In  addition,  the  data  sets  are  for  k1 . 


Figure  20:  Rotating  isotropic  decay  of  Shimomura.  Turbulent  kinetic  energy  versus  time. 
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In  Figures  21  and  22,  the  asterisks  represent  the  experimental  data  and  the  solid  lines 
represent  the  predictions  from  our  collision  model. 


Figure  21:  de  Bruyn  Kops  rotating  decaying  turbulence.  Turbulent  kinetic  energy  versus  time. 


t 


Figure  22:  Veeravalli's  decaying  kinetic  energy.  Kinetic  energy  versus  time. 
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In  Figure  23,  the  asterisks  represent  the  experimental  data  and  the  solid  lines  represent 
the  predictions  from  the  OEC  model. 
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Figure  23;  Rotating  isotropic  decay  of  Mansour,  Cambon  and  Speziale,  Turbulent  kinetic  energy 

versus  time,  a)  ReT=27.2  and  b)  Rer=67,t 
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Rapid  Distortion  Theory 


In  turbulent  shear  flows,  the  turbulence-to-mean-shear  time  scale  ratio  defined  as  SK/s 
varies  between  0  and  oo.  In  the  limiting  cases  when  the  ration  SK/e  is  exceptionally  large, 
the  evolution  of  the  turbulence  is  then  described  exactly  by  rapid-distortion  theory  or 
RDT.  Previous  work  compared  this  model  performance  to  that  of  a  standard  RDT  solver 
(i by  Chartrand).  This  time,  we  compare  our  model  performance  to  that  of  RDT  cases  of 
Matsumoto,  Blaisdell  and  Lee  &  Reynolds,  with  initial  conditions  summarized  in  Table7 
below.  Lee  &  Reynolds  experimented  three  cases:  axisymmetric  contraction  (AC), 
axisymmetric  expansion  (AE)  and  plane  strain  (PS).  Matsumoto’s  case  includes  two  DNS 
(high  and  low  Reynolds  numbers)  with  shear  (S)  deformation  while  Blaisdell  has  one 
elliptical  (E)  case. 


Lee  &  Reynolds 

(AQ (AE) (PS) 

Matsumoto 

(S) 

Blaisdell 

(E) 

e(m2/s3) 

0.018 

0.122 

0.25 

0.185 

1.79 

K(mV) 

1.0 

1.0 

1.0 

0.2 

1 

v(m2/s) 

10 

10 

10 

1.2e-2 

4.41e-2 

s  (s'1) 

1 

0.5 

1.0 

28.28 

3.0 

Rer 

5.59 

0.82 

0.4 

18.18 

12.75 

SK/e 

55.87 

4.08 

4 

30.6 

1.68 

Table  7:  Initial  conditions  of  Matsumoto  and  Lee&  Reynolds. 


Also,  included  in  Table  8  are  the  non -zero  mean  velocity  gradients  for  simple 
deformations. 


Axisymmetric 

contraction 

Axisymmetric 

expansion 

Plane 

Strain 

Shear 

5 

-2  5 

5 

0 

R22 

--S 

2 

5 

-S 

0 

-is 

2 

5 

0 

0 

*12 

0 

0 

0 

5 

5  =  (25*5,, O''2 

Ss 

2V35 

-25 

25 

Table  8:  Tensor  matrix  for  simple  deformations. 
The  graphs  below  summarize  the  results: 
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Figure  24:  Lee  &  Reynolds’  axisymmetric  contraction.  Tlie  dots  represent  the  DNS  and  the  lines 
represent  the  OEC  model  prediction,  SK/e=55,9, 


Figure  25:  Lee  &  Reynolds’  axisymmetric  expansion.  The  dots  represent  the  DNS  and  the  lines 
represent  the  OEC  model  prediction,  SK/e=9.08. 
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Figure  26:  Lee  &  Reynolds’  plane  strain,  The  dots  represent  the  DNS  and  the  lines  represent  the 

OEC  mode!  prediction.  SK/e=4. 


Matsumoto:  Rer=18.18 
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Figure  27:  Matsumoto’s  shear  deformation.  The  dots  represent  the  DNS  data  and  the  lines 
represent  the  OEC  model  prediction.  The  large  imposed  strain  (SK/e=30.6)  implies  RDT  is 

closely  approximated. 
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The  next  simulation  we  did  is  based  on  BlaisdelPs  DNS.  Here,  the  fact  that  both  the 
strain  ratio  and  the  turbulent  Reynolds  number  are  small  (respectively  1.68  and  12.75)  in 
addition  to  the  initial  random  field  justifies  the  RDT  approximation.  Furthermore,  we  ran 
four  simulations:  one  with  only  the  return-model  on,  a  second  one  with  just  the  rotation 
model  on,  a  third  one  with  both  return  and  rotation  models  on,  and  finally  the  RDT  case 
(return  and  rotation  models  both  turned  off).  Looking  at  the  graph  below,  we  were  able  to 
prove  that  both  return  and  rotation  have  no  effects  in  this  case.  The  results  are  shown 
below  in  Figure  28  and  29,  The  dots  represent  the  DNS  and  the  lines  show  the  model 
prediction. 


Figure  28:  BlaisdelPs  elliptical  flow  with  a)  return  model  on,  b)  rotation  model  on  and  c)  both 

return  and  rotation  models  on. 
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Figure  29:  Blaisdell's  elliptical  flow:  RDT 


Return-to-isotropy  Models 


For  anisotropic  cases,  a  term  to  model  the  return  to  isotropy  behavior  of  turbulent  flows 
was  introduced.  From  equation  (7a),  the  term  the  corresponds  to  the  return-to-isotropy 
model  for  the  Reynolds  stresses  is, 


(17) 


The  oriented-eddy  collision  model  includes  two  types  of  return  representations  H„  and 
k  -return. 


A 

Ry -return  model 


DtJ  was  modeled  in  the  following  ways: 


*  ft 

„  iLL,  " 


(18a) 

(18b) 

(18c) 

(18d) 

(18e) 


34 


The  first  two  equations  are  modeled  after  Rotta’s  Reynolds  Stress  Transport  (RST)  return 
models.  That  is,  both  equations  (18a)  and  (18b)  work  by  relaxing  each  individual 
Reynolds  stress  towards  an  isotropic  state  (e.g.  from  an  ellipse  to  a  sphere)  with  ( D,  ,B ) or 

without  (DtJA )  regard  to  the  other  eddies.  The  only  difference  between  the  two  equations 

is  that  one  uses  the  individual  kinetic  energy  of  each  eddy  ( K ),  while  the  second 
equation  uses  the  average  global  kinetic  energy  (■£);  thus  we  refer  to  D  A  as  Rotta-L 

(Local  Rotta)  and  DB  as  Rotta-G  (global  Rotta).  CR  is  a  tuned  constant  that  we 

determined  as  4.0  in  the  case  of  Rotta-L  and  2.5  for  Rotta-G.  Note  that  equations  (18c), 
(!8d)  and  (18e)  do  not  have  a  tunable  constant.  Those  were  generated  by  Perot  & 
Chartrand. 


/c-return  Models 

This  A-retum  model  is  part  of  the  orientation  equation  (7b)  and  corresponds  to: 


(19) 

The  term  ml  here  was  modeled  two  ways: 

(20a) 

with  and  *. -(*£*, A,)/(££k’) 

V=-Q!(3ffA'„ -sAk, 

(20b) 

withA^=x£(U/*2) 

The  first  equation  is  referred  to  as  the  “Ky  return  model”  while  the  second  one  is  the  “Ny 
return  model”.  Ny  depends  only  on  anisotropy  in  the  orientations  while  Ky  also  responds 
to  anisotropy  in  the  lengths  of  the  eddies.  CKl  and  CK2  are  tuning  constants  that  we 
determined  to  be  respectively  4.0  and  1.0.  Based  on  numerous  simulations,  we 
determined  Ky  to  be  the  best  performing  return  case  as  shown  below  in  Figure  30. 

At  first  it  seems  as  the  Kirreturn  model  performs  better  than  the  Ny,  However,  looking 
closely,  it  is  really  difficult  to  come  up  with  a  conclusion.  Ky  seems  to  work  best  on  the 
R„  terms  while  Ny  best  performs  on  the  non-diagonal  elements.  It  is  our  goal  to  further 
investigate  this  as  part  of  the  future  work. 
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Figure  30;  m*  and  mtH  model  comparisons 

As  previously  mentioned,  (CKi  ,CK2 )  are  tuning  constants  that  we  determined  to  be 
respectively  (4. 1 0)  for  Ny,  and  ( 1 ,4)  for  Ky. 


Shear/Strain  Flows 

In  this  section,  we  used  various  DNS  as  well  as  experimental  cases  to  test  our  mode! 
performance;  primary  in  the  Reynolds  stresses  analysis  of  shear  flows.  Tables  9  and  10 
below  provide  a  summary  with  the  values  of  the  constants  CR  andC*. : 


Matsumoto 

LePenven  A 

LePenven  B 

SK/e 

4.71 

0.43 

0.33 

Rct 

152 

612 

846 

(C, > CK ) 

(4,10) 

(4,10) 

(4,10) 

Strain  Tensor 

,r0  30  0^ 

0  0  0 

0  °J 

^5.48  0  0  ' 

0  1 .99  0 

v  0  0  -7.47  j 

^8.86  0  0  ' 

0  -2.36  0 

^  0  0  6.50 j 

Table  9:  Matsumoto  and  LePenven  summary  using  (Rotta-L,  KtJ),  which  is  (  DA  ,mA) 
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Hallback  -PS 

Rej 

11 

»  Of  ) 

0,4) 

SK/e 

9 

3 

1 

Strain 

Tensor 

( 4.36  0  0^ 

0  -4.36  0 

,  0  0  0, 

0.46  0  0  N 

0  0  0 
v  0  0  -1.46 ; 

r0.49  0  0  ' 

0  0  0 
k  0  0  -0.49; 

Table  10:  Hallback’s  summary  using  (Rotta-L,  K;j)  for  Plane  Strain 


Numerical  Results:  return-to-isotropy  and  shear/strain  deformation 

To  illustrate  the  retum-to-isotropy  above,  three  different  cases  were  used:  LePenven, 
Matsumoto  and  Hallback.  Only  cases  that  are  highly  dependent  on  return-to-isotropy 
were  used  to  determine  the  values  of  our  constants  as  well  as  to  validate- the  models.  All 
initial  conditions  are  shown  above  in  Tables  9  and  10.  The  results  are  shown  in  Figures 
31  and  32. 


figure  3 1 :  Le  Penven  -  case  A.  a)  Reynolds  stresses  and  b)  Kinetic  energy 


Figure  32:  Le  Penven  -  case  B.  a)  Reynolds  stresses  and  b)  Kinetic  energy 
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Correspondingly,  the  oriented-eddy  collision  prediction  was  compared  to  the 
homogeneous  shear  and  strain  flows:  Matsumoto  (Figure  33),  and  Hallback  PS  (Figure 
34): 


Figure  33:  Matsumoto’s  shear  deformation.  The  dots  represent  his  DNS  and  the  lines  represent 
our  model  prediction,  a)  Reynolds  stresses  and  b)  Dissipation  and  kinetic  energy 
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Figure  34:  Hallback  -  Plane  Strain  a)  S=3  and  b)  S=1 


39 


Diffusion 

In  equations  (7a)  and  (7b)  the  final  terms,  v|(v  +  t/r)V.fltf  J  models  the  diffusive  action  of 
the  Reynolds  stresses  while  V((v  +  v7.)VA,)  accounts  for  the  diffusive  action  of  the 
orientation  vectors  .  In  one-dimension,  V^(v  +  vT )Vtf1(  J  corresponds  to 


d  *  ,  dR 

— (y+y.) — - 

dyX  dy 


(21) 


v  is  the  fluid  viscosity  while  v,  corresponds  to  the  eddy  viscosity.  We  defined  local  and 
global  eddy  viscosities.  As  mentioned  before,  “local’'  implies  that  all  calculations  are 
done  locally.  In  this  case,  the  model  uses  a  local  v,  that  is  defined  as 


»  I  A 

with  K  =  —  Rn  and  CL=  1 .  Regarding  the  global  eddy  viscosities,  two  equations  that 

are  referred  to  as  globabll  and  global2  (V/'1  and  v/'2 )  are  currently  being  evaluated.  The 

global  1  and  global2  are  similar,  with  the  only  difference  that  in  the  first  case,  the  local 
kinetic  energy  and  dissipation  are  being  summed  before  dividing,  whereas  in  the  second 
case,  the  summation  is  done  after  the  division.  These  concepts  are  illustrated  below 
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(22b) 


(22c) 


where  CGl  -  CG2  =  1  and  E  implies  summation  over  the  orientations. 


As  expected  tor  isotropic  flows  (the  orientations  vectors  all  have  then  same  length),  both 
global  eddy  viscosity  formulas  {vf]  and  v,0*)  performed  equally  as  shown  in  Figure  35 
below.  The  comparison  was  done  using  the  DNS  of  Carati: 
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Carati:  vj-  comparison 
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Figure  35:  Eddy  viscosity  comparison  for  both  global  equations 


After  implementing  the  diffusion  in  the  code  along  with  all  three  variants  of  the  eddy 
viscosity  various  simulations  were  conducted  in  order  to  determine  the  efficacy  of  the 
eddy  collision  model.  It  is  important  to  mention  that  the  kinetic  energy  decay  is  no 
longer  homogeneous  (as  previously)  but  instead  is  also  spatially  dependent. "in  the 
diffusion  case,  at  one  fixed  time  t,  we  are  looking  at  both  the  kinetic  energy  and 
dissipation  at  different  locations  (y).  The  first  step  in  the  analysis  is  to  determine  which 
eddy  viscosity  equation  best  models  the  diffusion  process.  Starting  with  i/(  01  (22a), 

various  simulations  were  conducted  as  part  of  the  evaluation  process.  The  first  simulation 
was  run  against  that  of  Chasnov  and  shows  the  diffusion  process  at  different  times  t. 
Chasnov’s  flow  is  inhomogeneous  with  the  following  characteristics:  shearless, 
irrotational  and  isotropic  with  periodic  boundary  conditions.  Note  that  in  order  to  reduce 
the  time  step,  we  interpolated  the  original  data  as  represented  by  the  solid  blue  line. 

Next  in  Figure  36,  we  looked  at  the  diffusion  evolution  at  times  t=  1.375,  4.125  and  9.625 
seconds.  The  asterisks  represent  the  data  and  the  matching  solid  blue  lines  correspond  to 
our  simulations. 
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Figure  36:  Kinetic  Energy  versus  position  at  different  times  t.  The  matching  blue  lines  correspond 
to  the  OEC  simulations.  a)li near-linear  plot  and  b)log-l  inear  plot. 
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The  second  diffusion  simulation  matched  that  of  Gilbert.  Gilbert  assumes  a  shearless, 
irrotational  and  homogeneous  flow.  In  addition,  the  flow  has  some  levels  of  anisotropy. 
The  stars  represent  data  from  Gilbert  at  times  t=0,  0.0292,  0.0402,  0.0764,  0.0884, 
0.1 154,  0.1274,  0.1634  and  0.2024  seconds.  The  matching  solid  lines  correspond  to  the 
simulations. 
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Figure  37:  RM  (kinetic  energy  component)  versus  position.  The  stars  represent  data  from  Gilbert 
at  times  t=0,  0.0292,  0.0402,  0.0764,  0.0884,  0.1 154,  0.1274,  0.1634  and  0.2024  seconds.  The 
matching  solid  lines  correspond  to  the  OEC  simulations. 


The  final  set  of  data  that  was  looked  at  is  a  more  recent  one  (2003)  and  was  published  by 
Carati.  Carati's  data  is  unique  in  a  sense  that  we  have  access  to  both  the  kinetic  energy 
and  the  dissipation  rate.  Here  although  not  ideal,  we  used  zero  boundary  conditions 
compared  to  periodic  conditions  in  the  two  cases  above  (Chasnov,  Gilbert).  For  reasons 
that  remain  unclear  at  this  time,  the  OEC  isotropic  simulations  decay  a  little  faster  than 
expected.  The  results  obtained  are  shown  below  in  Figures  38  and  39: 
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Figure  38;  Kinetic  Energy  versus  position  at  different  times  t.  The  stars  represent  data  from  Carati 
at  times  t=0,  0.071  and  0. 191  seconds.  The  matching  solid  lines  correspond  to  our  simulations. 


Figure  39:  Dissipation  versus  position  at  different  times  t.  The  stars  represent  data  from  Carati  at 
times  1=0,  0.071  and  0. 191  seconds.  The  matching  solid  lines  correspond  to  our  simulations. 
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Conclusions  &  Future  Work 

This  project  has  allowed  us  to  demonstrate  that  oriented-eddy  collisional  (OEC)  models 
are  an  interesting,  accurate,  and  viable  approach  to  turbulence  modeling.  We  have 
demonstrated  that: 

•  Models  exist  in  the  regime  between  LES  and  RANS  that  have  very  attractive  cost 
and  accuracy  attributes  for  current  day  design. 

•  It  is  possible  to  increase  the  physics  in  turbulence  models  and  reduce  the  number 
of  tuned  constants,  while  still  having  a  cost  effective  model  that  can  run  on  a  PC. 

•  The  structure  (orientation)  of  turbulence  is  just  as  import  as  the  magnitude  of  the 
fluctuations.  Models  which  represent  structure  have  huge  advantages  in  capturing 
the  turbulence  physics. 

•  The  model  can  be  interpreted  as  a  model  for  the  evolution  of  the  two-point 
correlation.  Critical  to  this  model  —  is  decomposing  the  two-point  correlation 
into  self-similar  ‘modes’. 

As  with  any  turbulence  model,  a  great  deal  of  work  remains  to  validate  this  model.  In 
this  project  we  have  clearly  demonstrated  that  the  approach  is  extensible  and  can 

accurately  predict  a  wide  variety  of  quite  different  but  fundamental  turbulent  flow 
situations. 

Future  work  will  complete  the  modeling  of  wall  effects.  In  addition,  we  expect  this 
model  to  predict  transition  very  well,  and  this  will  be  demonstrated,  Finally,  this  model 
will  be  implemented  in  a  3D,  unstructured,  parallel,  Navier-Stokes  code  so  that  more 
complex  and  practical  flow  situations  can  be  tested. 
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